This vignette shows some examples on how to plot pretty maps in R with the latest packages. Currently, it functions as a loose collection of open-source and third-party map material combined with various raster and shapefiles. I would already like to mention that the most beautiful map was completely taken from here: https://timogrossenbacher.ch/2019/04/bivariate-maps-with-ggplot2-and-sf/ Generally, all code chunks in this vignette need to run in tsa with the current R-version 3.5.2 (2020-11-20).
Quick Note about the setup: I access Tsa@CSCS via VSCode and remote-ssh. Then I run an interactive R sessions while working in the vignette. This required quite a few steps to set it up properly. I plan to write a wiki page about vscode-R @ cscs at some point in the future. Feel free to ask me anything in the meantime. In the chunk below I load a whole lot of packages, if the reader is only interested in a portion of the plot, feel free to only load a selection of packages as well.
A typical output format in weather and climate science are netcdf files, hence we import one of this type right here, showing the median pollen concentrations modeled by Cosmo-1E. The code can of course be adjusted to plot any netcdf file.
Different packages will require different formats. The most recent development when it comes to polygon data in R is the package sf, a highly flexible framework not yet supported by all other packages. Hence we prepare the data in different shapes and forms. The output of many weather models is a gridded field. First, we display such grids. In a second step we extract data from model output for specific locations (Swiss municipalities) to create even nicer plots specifically for the Swiss Domain.
con <- nc_open(nc_path)
layer <- ncvar_get(con, "POAC")[, , 80]
x <- ncvar_get(con, "lon_1")
y <- ncvar_get(con, "lat_1")
nc_close(con)
r <- raster(layer,
ymn = 42.67, ymx = 49.52,
xmn = 0.16, xmx = 16.75, crs = "+proj=longlat"
)
layer_latlong <- tibble(
x = c(x),
y = c(y),
layer = c(layer)
) %>%
# lons west of 0 deg. need to be negative, hence make them negative here
mutate(x = if_else(x > 180, x - 360, x))
# Convert the banana-shape into an equally spaced grid
layer_raster <- rasterize(cbind(c(x), c(y)),
r, layer_latlong$layer,
fun = mean
)
# Plots can become heavy - adjust second agrument accordingly
layer_raster_vcoarse <- aggregate(layer_raster, 20)
layer_raster_coarse <- aggregate(layer_raster, 6)
# Polygons are usually faster to plot than raster images
layer_poly <- rasterToPolygons(layer_raster_coarse)
layer_poly_leaflet <- rasterToPolygons(layer_raster_vcoarse)
# This is the oldschool way before sf was developped
layer_poly@data$id <- seq_len(nrow(layer_poly@data))
poly_fort <- tidy(layer_poly, data = layer_poly@data)
# join data
poly_fort_mer <- merge(poly_fort, layer_poly@data,
by.x = "id", by.y = "id"
)
# This is the modern way and works well with ggmap
layer_poly_sf <- st_as_sf(layer_poly)
We first define a general map-them (from Timo Grossenbacher) to obtain a similar look and feel for all plots:
The following creates an interactive map (html widget) where the background can be chosen from a large variety of map providers. If the polygon layer is fine or dense, plotting can take a long time or even crash. It’s usually better to first aggregate the polygons.
pal <- colorNumeric("YlOrRd", layer_poly_leaflet$layer)
mymap <- "https://stamen-tiles-{s}.a.ssl.fastly.net/terrain/{z}/{x}/{y}{r}.png"
leaflet_map <- leaflet(layer_poly_leaflet) %>%
# More maps here: https://leaflet-extras.github.io/leaflet-providers/preview/
addTiles(urlTemplate = mymap) %>%
# Cosmo 1E Domain with extended boundaries
fitBounds(lat1 = 42.60, lat2 = 49.60, lng1 = 0.10, lng2 = 16.80) %>%
addPolygons(
weight = 0,
popup = as.character(round(values(layer_raster_vcoarse), 2)),
smoothFactor = 0.5, fillColor = ~ pal(layer), fillOpacity = 0.6
) %>%
addLegend(pal = pal, values = ~ layer_poly_leaflet$layer, title = "Poaceae")
leaflet_map
The following plots are fully open source.
upper_left <- c(51, 0)
lower_right <- c(40, 17.5)
# osm and other map providers available
map <- openmap(upper_left, lower_right, type = "stamen-terrain")
# The raster map is in a mercator projection and must be transformed
map_proj <- openproj(map)
raster_map <- autoplot(map_proj) +
my_maptheme() +
geom_polygon(
data = poly_fort_mer,
aes(x = long, y = lat, group = group, fill = layer),
alpha = 0.7,
size = 0
) +
scale_fill_gradientn("Poaceae\n[m^-3]", colours = rev(heat.colors(200))) +
labs(
x = "Lon",
y = "Lat",
title = "Grass Pollen in Cosmo-1 Domain",
subtitle = "Hourly Average Concentration on the 1st of July 2020 at Midnight
"
)
raster_map
The next example requires a Google account to retrieve data from the Google API. This package is quite powerful and works well with the modern sf datatype.
# ?register_google
# Enter your key here - I will change mine soon ;-)
register_google(key = "AIzaSyC-8Tyjbj1zc5fKdKZdbWaJk1H8CXdUIB8")
centroid <- c(45, 8)
lat_zoom <- c(45, 48)
long_zoom <- c(4, 11)
ggmap_heat <- get_map(c(lon = centroid[2], lat = centroid[1]),
zoom = 6, maptype = "terrain", color = "color", scale = "auto"
) %>%
ggmap() +
geom_sf(aes(fill = layer, alpha = layer),
data = layer_poly_sf, inherit.aes = FALSE, lwd = 0
) +
my_maptheme() +
theme(panel.grid.major = element_line(color = "white")) +
scale_fill_gradientn("Grass\nPollen", colors = rev(heat.colors(200))) +
scale_alpha_continuous(range = c(0.1, 0.9)) +
scale_x_continuous(limits = long_zoom, expand = c(0, 0)) +
scale_y_continuous(limits = lat_zoom, expand = c(0, 0)) +
guides(alpha = "none") +
labs(
x = "Lon",
y = "Lat",
title = "Grass Pollen in Switzerland",
subtitle = "Hourly Average Concentration on the 1st of July 2020 at Midnight
"
)
ggmap_heat